home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / docs / apps / circuits / spice2g6.z / spice2g6 / spice / Fortran / moseq2.f < prev    next >
Encoding:
Text File  |  1989-02-03  |  12.7 KB  |  447 lines
  1.       subroutine moseq2(vds,vbs,vgs,gm,gds,gmbs,
  2.      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
  3.       implicit double precision (a-h,o-z)
  4. c
  5. c     this routine evaluates the drain current, its derivatives and
  6. c     the charges associated with the gate, channel and bulk
  7. c     for mosfets
  8. c
  9. c spice version 2g.6  sccsid=mosarg 3/15/83
  10.       common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld,
  11.      1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof,
  12.      2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev
  13. c spice version 2g.6  sccsid=status 3/15/83
  14.       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
  15.      1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
  16.      2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
  17. c spice version 2g.6  sccsid=knstnt 3/15/83
  18.       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
  19.      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
  20.      2   pivtol,pivrel
  21. c
  22.       dimension a4(4),b4(4),x4(8),poly4(8),sig1(4),sig2(4)
  23.       data sig1 / 1.0d0, -1.0d0, 1.0d0, -1.0d0/,
  24.      1     sig2 / 1.0d0,  1.0d0,-1.0d0, -1.0d0/
  25. c
  26. c     icharg=1 causes charges to be computed
  27. c     icharg=0 bypasses the computation of charges
  28. c
  29.       icharg=1
  30.       if (mode.ne.1.and.xqco.le.0.5d0) go to 100
  31.       icharg=0
  32.       if (xqco.gt.0.5d0) go to 100
  33.       if (modedc.eq.2.and.nosolv.ne.0) icharg=1
  34.       if (initf.eq.4) icharg=1
  35. c
  36. c  compute some useful quantities
  37. c
  38.   100 if (vbs.gt.0.0d0) go to 110
  39.       sarg=dsqrt(phi-vbs)
  40.       tsarg=sarg+sarg
  41.       dsrgdb=-0.5d0/sarg
  42.       d2sdb2=+0.5d0*dsrgdb/(phi-vbs)
  43.       go to 120
  44.   110 sphi=dsqrt(phi)
  45.       sphi3=phi*sphi
  46.       sarg=sphi/(1.0d0+0.5d0*vbs/phi)
  47.       tsarg=sarg+sarg
  48.       dsrgdb=-0.5d0*sarg*sarg/sphi3
  49.       d2sdb2=-dsrgdb*sarg/sphi3
  50.   120 if ((vds-vbs).lt.0.0d0) go to 130
  51.       barg=dsqrt(phi+vds-vbs)
  52.       dbrgdb=-0.5d0/barg
  53.       d2bdb2=+0.5d0*dbrgdb/(phi+vds-vbs)
  54.       go to 200
  55.   130 barg=sphi/(1.0d0+0.5d0*(vbs-vds)/phi)
  56.       dbrgdb=-0.5d0*barg*barg/sphi3
  57.       d2bdb2=-dbrgdb*barg/sphi3
  58. c
  59. c  calculate threshold voltage (von)
  60. c     narrow-channel effect
  61. c
  62.   200 factor=0.125d0*fnarrw*twopi*epssil/cox*xl
  63.       eta=1.0d0+factor
  64.       vbin=vbi+factor*(phi-vbs)
  65.       if (gamma.le.0.0d0) go to 215
  66.       if (xnsub.le.0.0d0) go to 215
  67.       xwd=xd*barg
  68.       xws=xd*sarg
  69. c
  70. c     short-channel effect with vds .ne. 0.0d0
  71. c
  72.       argss=0.0d0
  73.       argsd=0.0d0
  74.       dbargs=0.0d0
  75.       dbargd=0.0d0
  76.       dgdvds=0.0d0
  77.       dgddb2=0.0d0
  78.       if (xj.le.0.0d0) go to 205
  79.       argxs=1.0d0+2.0d0*xws/xj
  80.       args=dsqrt(argxs)
  81.       argss=0.5d0*xj/xl*(args-1.0d0)
  82.       argxd=1.0d0+2.0d0*xwd/xj
  83.       argd=dsqrt(argxd)
  84.       argsd=0.5d0*xj/xl*(argd-1.0d0)
  85.   205 gamasd=gamma*(1.0d0-argss-argsd)
  86.       gamass=gamma*(1.0d0-2.0d0*argss)
  87.       dbxwd=xd*dbrgdb
  88.       dbxws=xd*dsrgdb
  89.       if (xj.le.0.0d0) go to 210
  90.       dbargs=0.5d0/xl*dbxws/args
  91.       dbargd=0.5d0/xl*dbxwd/argd
  92.       dasdb2=-xd*( d2sdb2+dsrgdb*dsrgdb*xd/(xj*argxs) )/(xl*args)
  93.       daddb2=-xd*( d2bdb2+dbrgdb*dbrgdb*xd/(xj*argxd) )/(xl*argd)
  94.       dgddb2=-0.5d0*gamma*(dasdb2+daddb2)
  95.   210 dgddvb=-gamma*(dbargs+dbargd)
  96.       dgsdvb=-2.0d0*gamma*dbargs
  97.       if (xj.le.0.0d0) go to 220
  98.       ddxwd=-dbxwd
  99.       dgdvds=-gamma*0.5d0/xl*ddxwd/argd
  100.       go to 220
  101.   215 gamasd=gamma
  102.       gamass=gamma
  103.       gammad=gamma
  104.       dgddvb=0.0d0
  105.       dgsdvb=0.0d0
  106.       dgdvds=0.0d0
  107.       dgddb2=0.0d0
  108.   220 von=vbin+gamasd*sarg
  109. c     write(iofile,221) von,vbin,vbi,gamasd,argss,argsd,xj
  110.   221 format ('0msg1:'/1p7d10.2)
  111.       vth=von
  112.       vdsat=0.0d0
  113.   225 if (xnfs.eq.0.0d0.or.cox.eq.0.0d0) go to 230
  114.       cfs=charge*xnfs
  115.       cdonco=-(gamasd*dsrgdb+dgddvb*sarg)+factor
  116.       xn=1.0d0+cfs/cox*xw*xl+cdonco
  117.       von=von+vt*xn
  118. c     write (iofile,226) von,cdonco,xn,cfs,xd
  119.   226 format(' msg2:'/1p6d10.2)
  120.       argg=1.0d0/(vt*xn)
  121.       vgst=vgs-von
  122.       go to 300
  123.   230 vgst=vgs-von
  124.       if (vgs.gt.von) go to 300
  125. c
  126. c  cutoff region
  127. c
  128.       gds=0.0d0
  129.       go to 1050
  130. c
  131. c  compute some more useful quantities
  132. c
  133.   300 sarg3=sarg*sarg*sarg
  134.       sbiarg=dsqrt(phib)
  135.       gammad=gamasd
  136.       dgdvbs=dgddvb
  137.       body=barg*barg*barg-sarg3
  138.       gdbdv=2.0d0*gammad*(barg*barg*dbrgdb-sarg*sarg*dsrgdb)
  139.       dodvbs=-factor+dgdvbs*sarg+gammad*dsrgdb
  140.       if (xnfs.eq.0.0d0) go to 400
  141.       if (cox.eq.0.0d0) go to 410
  142.       dxndvb=2.0d0*dgdvbs*dsrgdb+gammad*d2sdb2+dgddb2*sarg
  143.       dodvbs=dodvbs+vt*dxndvb
  144.       dxndvd=dgdvds*dsrgdb
  145.       dodvds=dgdvds*sarg+vt*dxndvd
  146. c
  147. c  evaluate effective mobility and its derivatives
  148. c
  149.   400 if (cox.le.0.0d0) go to 410
  150.       udenom=vgst
  151.       if (udenom.le.vbp) go to 410
  152.       ufact=dexp(uexp*dlog(vbp/udenom))
  153.       ueff=uo*ufact
  154.       dudvgs=-ufact*uexp/udenom
  155.       dudvds=0.0d0
  156.       dudvbs=uexp*ufact*dodvbs/vgst
  157.       go to 500
  158.   410 ufact=1.0d0
  159.       ueff=uo
  160.       dudvgs=0.0d0
  161.       dudvds=0.0d0
  162.       dudvbs=0.0d0
  163. c
  164. c     evaluate saturation voltage and its derivatives according to
  165. c     grove-frohman equation
  166. c
  167.   500 vgsx=vgs
  168.       gammad=gamasd/eta
  169.       dgdvbs=dgddvb
  170.       if (xnfs.ne.0.0d0.and.cox.ne.0.0d0)
  171.      1   vgsx=dmax1(vgs,von)
  172.   505 if (gammad.le.0.0d0) go to 535
  173.       gammd2=gammad*gammad
  174.       argv=(vgsx-vbin)/eta+phi-vbs
  175.       if (argv.le.0.0d0) go to 540
  176.       arg=dsqrt(1.0d0+4.0d0*argv/gammd2)
  177.       vdsat=(vgsx-vbin)/eta+gammd2*(1.0d0-arg)/2.0d0
  178.       vdsat=dmax1(vdsat,0.0d0)
  179.   510 if (icharg.eq.0) go to 530
  180.       arg1=gammd2/(eta*eta)
  181.       arg2=vds-0.5d0*arg1
  182.       argsq=(arg2+0.5d0*arg1+phi-vbs)*arg1
  183.       if (argsq.ge.0.0d0) go to 515
  184.       vpof=vth
  185.       go to 520
  186.   515 vpof=vbin+eta*(arg2+0.5d0*arg1+dsqrt(argsq))
  187.   520 argv1=(vpof-vbin)/eta+phi-vbs
  188.       if (argv1.gt.0.0d0) go to 525
  189.       vdsat1=0.0d0
  190.       go to 530
  191.   525 arg1=dsqrt(1.0d0+4.0d0*argv1/gammd2)
  192.       vdsat1=(vpof-vbin)/eta+gammd2*(1.0d0-arg1)/2.0d0
  193.       vdsat1=dmax1(vdsat1,0.0d0)
  194.   530 dsdvgs=(1.0d0-1.0d0/arg)/eta
  195.       dsdvbs=(gammad*(1.0d0-arg)+2.0d0*argv/(gammad*arg))/eta*dgdvbs+
  196.      1       1.0d0/arg+factor*dsdvgs
  197.       go to 545
  198.   535 vdsat=dmax1((vgsx-vbin)/eta,0.0d0)
  199.       vpof=dmax1((eta*vds+vbin),0.0d0)
  200.       vdsat1=dmax1((vpof-vbin)/eta,0.0d0)
  201.       dsdvgs=1.0d0
  202.       dsdvbs=0.0d0
  203.       go to 545
  204.   540 vdsat=0.0d0
  205.       vpof=vth
  206.       vdsat1=0.0d0
  207.       dsdvgs=0.0d0
  208.       dsdvbs=0.0d0
  209. c
  210. c     store vdsat as above in vpof (pinch-off)
  211. c
  212.   545 if (vmax.le.0.0d0) go to 600
  213. c
  214. c     evaluate saturation voltage and its derivatives according to
  215. c     baum's theory of scattering velocity saturation
  216. c
  217.       gammd2=gammad*gammad
  218.       v1=(vgsx-vbin)/eta+phi-vbs
  219.       v2=phi-vbs
  220.       xv=vmax*xl/ueff
  221.       a1=gammad/0.75d0
  222.       b1=-2.0d0*(v1+xv)
  223.       c1=-2.0d0*gammad*xv
  224.       d1=2.0d0*v1*(v2+xv)-v2*v2-4.0d0/3.0d0*gammad*sarg3
  225.       a=-b1
  226.       b=a1*c1-4.0d0*d1
  227.       c=-d1*(a1*a1-4.0d0*b1)-c1*c1
  228.       r=-a*a/3.0d0+b
  229.       s=2.0d0*a*a*a/27.0d0-a*b/3.0d0+c
  230.       r3=r*r*r
  231.       s2=s*s
  232.       p=s2/4.0d0+r3/27.0d0
  233.       p0=dabs(p)
  234.       p2=dsqrt(p0)
  235.       if (p.ge.0.0d0) go to 550
  236.       ro=dsqrt(s2/4.0d0+p0)
  237.       ro=dlog(ro)/3.0d0
  238.       ro=dexp(ro)
  239.       fi=datan(-2.0d0*p2/s)
  240.       y3=2.0d0*ro*dcos(fi/3.0d0)-a/3.0d0
  241.       go to 560
  242.   550 p3=dexp(dlog(dabs(-s/2.0d0+p2))/3.0d0)
  243.       p4=dexp(dlog(dabs(-s/2.0d0-p2))/3.0d0)
  244.       y3=p3+p4-a/3.0d0
  245.   560 iknt=0
  246.       a3=dsqrt(a1*a1/4.0d0-b1+y3)
  247.       b3=dsqrt(y3*y3/4.0d0-d1)
  248.       do 570 i=1,4
  249.       a4(i)=a1/2.0d0+sig1(i)*a3
  250.       b4(i)=y3/2.0d0+sig2(i)*b3
  251.       delta4=a4(i)*a4(i)/4.0d0-b4(i)
  252.       if (delta4.lt.0.0d0) go to 570
  253.       iknt=iknt+1
  254.       x4(iknt)=-a4(i)/2.0d0+dsqrt(delta4)
  255.       iknt=iknt+1
  256.       x4(iknt)=-a4(i)/2.0d0-dsqrt(delta4)
  257.   570 continue
  258.       jknt=0
  259.       do 580 j=1,iknt
  260.       if (x4(j).le.0.0d0) go to 580
  261.       poly4(j)=x4(j)*x4(j)*x4(j)*x4(j)+a1*x4(j)*x4(j)*x4(j)
  262.       poly4(j)=poly4(j)+b1*x4(j)*x4(j)+c1*x4(j)+d1
  263.       if (dabs(poly4(j)).gt.1.0d-6) go to 580
  264.       jknt=jknt+1
  265.       if (jknt.gt.1) go to 575
  266.       xvalid=x4(j)
  267.   575 if (x4(j).gt.xvalid) go to 580
  268.       xvalid=x4(j)
  269.   580 continue
  270.       if (jknt.gt.0) go to 590
  271.       ivmflg=ivmflg+1
  272.       go to 600
  273.   590 vdsat=xvalid*xvalid+vbs-phi
  274. c
  275. c  evaluate effective channel length and its derivatives
  276. c
  277.   600 if (vds.eq.0.0d0) go to 610
  278.       gammad=gamasd
  279.       if ((vbs-vdsat).gt.0.0d0) go to 601
  280.       bsarg=dsqrt(vdsat-vbs+phi)
  281.       dbsrdb=-0.5d0/bsarg
  282.       go to 602
  283.   601 bsarg=sphi/(1.0d0+0.5d0*(vbs-vdsat)/phi)
  284.       dbsrdb=-0.5d0*bsarg*bsarg/sphi3
  285.   602 bodys=bsarg*bsarg*bsarg-sarg3
  286.       gdbdvs=2.0d0*gammad*(bsarg*bsarg*dbsrdb-sarg*sarg*dsrgdb)
  287.       if (vmax.gt.0.0d0) go to 603
  288.       if (xnsub.eq.0.0d0) go to 610
  289.       if (xlamda.gt.0.0d0) go to 610
  290.       argv=(vds-vdsat)/4.0d0
  291.       sargv=dsqrt(1.0d0+argv*argv)
  292.       arg=dsqrt(argv+sargv)
  293.       xlfact=xd/(xl*vds)
  294.       xlamda=xlfact*arg
  295.       dldsat=vds*xlfact*arg/(8.0d0*sargv)
  296.       go to 605
  297.   603 argv=(vgsx-vbin)/eta-vdsat
  298.       xdv=xd/dsqrt(xneff)
  299.       xlv=vmax*xdv/(2.0d0*ueff)
  300.       vqchan=argv-gammad*bsarg
  301.       dqdsat=-1.0d0+gammad*dbsrdb
  302.       vl=vmax*xl
  303.       dfunds=vl*dqdsat-ueff*vqchan
  304.       dfundg=(vl-ueff*vdsat)/eta
  305.       dfundb=-vl*(1.0d0+dqdsat-factor/eta)+
  306.      1        ueff*(gdbdvs-dgdvbs*bodys/1.5d0)/eta
  307.       dsdvgs=-dfundg/dfunds
  308.       dsdvbs=-dfundb/dfunds
  309.       if (xnsub.eq.0.0d0) go to 610
  310.       if (xlamda.gt.0.0d0) go to 610
  311.       argv=dmax1(vds-vdsat,0.0d0)
  312.       xls=dsqrt(xlv*xlv+argv)
  313.       dldsat=xdv/(2.0d0*xls)
  314.       xlfact=xdv/(xl*vds)
  315.       xlamda=xlfact*(xls-xlv)
  316.       dldsat=dldsat/xl
  317.   605 dldvgs=dldsat*dsdvgs
  318.       dldvds=-xlamda+dldsat
  319.       dldvbs=dldsat*dsdvbs
  320.       go to 620
  321.   610 dldvgs=0.0d0
  322.       dldvds=0.0d0
  323.       dldvbs=0.0d0
  324. c
  325. c     limit channel shortening at punch-through
  326. c
  327.   620 xwb=xd*sbiarg
  328.       xld=xl-xwb
  329.       clfact=1.0d0-xlamda*vds
  330.       dldvds=-xlamda-dldvds
  331.       xleff=xl*clfact
  332.       deltal=xlamda*vds*xl
  333.       if (xnsub.eq.0.0d0) xwb=0.25d-6
  334.       if (xleff.ge.xwb) go to 700
  335.       xleff=xwb/(1.0d0+(deltal-xld)/xwb)
  336.       clfact=xleff/xl
  337.       dfact=xleff*xleff/(xwb*xwb)
  338.       dldvgs=dfact*dldvgs
  339.       dldvds=dfact*dldvds
  340.       dldvbs=dfact*dldvbs
  341. c
  342. c  evaluate effective beta (effective kp)
  343. c
  344.   700 beta1=beta*ufact/clfact
  345. c
  346. c  test for mode of operation and branch appropriately
  347. c
  348.       gammad=gamasd
  349.       dgdvbs=dgddvb
  350.       if (vds.gt.1.0d-8) go to 730
  351.       if (vgs.gt.von) go to 720
  352.       if ((xnfs.ne.0.0d0).and.(cox.ne.0.0d0)) go to 710
  353.       gds=0.0d0
  354.       go to 1050
  355. c
  356.   710 gds=beta1*(von-vbin-gammad*sarg)*dexp(argg*(vgs-von))
  357.       go to 1050
  358. c
  359. c
  360.   720 gds=beta1*(vgs-vbin-gammad*sarg)
  361.       go to 1050
  362. c
  363.   730 if (vgs.gt.von) go to 900
  364. c
  365. c  subthreshold region
  366. c
  367.       if (vdsat.gt.0.0d0) go to 830
  368.       gds=0.0d0
  369.       if (vgs.gt.vth) go to 1020
  370.       go to 1050
  371.   830 vdson=dmin1(vdsat,vds)
  372.       if (vds.le.vdsat) go to 850
  373.       barg=bsarg
  374.       dbrgdb=dbsrdb
  375.       body=bodys
  376.       gdbdv=gdbdvs
  377.   850 cdson=beta1*((von-vbin-eta*vdson*0.5d0)*vdson-gammad*body/1.5d0)
  378.       didvds=beta1*(von-vbin-eta*vdson-gammad*barg)
  379.       gdson=-cdson*dldvds/clfact-beta1*dgdvds*body/1.5d0
  380.       if (vds.lt.vdsat) gdson=gdson+didvds
  381.       gbson=-cdson*dldvbs/clfact
  382.      1      +beta1*(dodvbs*vdson+factor*vdson-dgdvbs*body/1.5d0-gdbdv)
  383.       if (vds.gt.vdsat) gbson=gbson+didvds*dsdvbs
  384.       expg=dexp(argg*(vgs-von))
  385.       cdrain=cdson*expg
  386.       gmw=cdrain*argg
  387.       gm=gmw
  388.       if (vds.gt.vdsat) gm=gmw+didvds*dsdvgs*expg
  389.       gds=gdson*expg-gm*dodvds-gmw*(vgs-von)*dxndvd/xn
  390.       gmbs=gbson*expg-gm*dodvbs-gmw*(vgs-von)*dxndvb/xn
  391.       go to 1020
  392. c
  393. c
  394.   900 if (vds.gt.vdsat) go to 1000
  395. c
  396. c  linear region
  397. c
  398.       cdrain=beta1*((vgs-vbin-eta*vds/2.0d0)*vds-gammad*body/1.5d0)
  399.       arg=cdrain*(dudvgs/ufact-dldvgs/clfact)
  400.       gm=arg+beta1*vds
  401.       arg=cdrain*(dudvds/ufact-dldvds/clfact)
  402.       gds=arg+beta1*(vgs-vbin-eta*vds-
  403.      1   gammad*barg-dgdvds*body/1.5d0)
  404.       arg=cdrain*(dudvbs/ufact-dldvbs/clfact)
  405.       gmbs=arg-beta1*(gdbdv+dgdvbs*body/1.5d0-factor*vds)
  406.       go to 1020
  407. c
  408. c  saturation region
  409. c
  410.  1000 cdrain=beta1*((vgs-vbin-eta*vdsat/2.0d0)*vdsat-gammad*bodys/1.5d0)
  411.       arg=cdrain*(dudvgs/ufact-dldvgs/clfact)
  412.       gm=arg+beta1*vdsat+
  413.      1   beta1*(vgs-vbin-eta*vdsat-gammad*bsarg)*dsdvgs
  414.       gds=-cdrain*dldvds/clfact-beta1*dgdvds*bodys/1.5d0
  415.       arg=cdrain*(dudvbs/ufact-dldvbs/clfact)
  416.       gmbs=arg-beta1*(gdbdvs+dgdvbs*bodys/1.5d0-factor*vdsat)+
  417.      1     beta1*(vgs-vbin-eta*vdsat-gammad*bsarg)*dsdvbs
  418. c
  419. c     compute charges for "on" region
  420. c
  421.  1020 if (icharg.eq.0) go to 1500
  422.       if (vgs.le.vth) go to 1070
  423.       call mqspof(vds,vbs,vgs,vpof,vdsat1,vth,vbin,gamasd,
  424.      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
  425.       go to 2000
  426. c
  427. c  finish special cases
  428. c
  429.  1050 cdrain=0.0d0
  430.       gm=0.0d0
  431.       gmbs=0.0d0
  432.  1070 xqc=xqco
  433.       if (icharg.eq.0) go to 1500
  434.       call mosq2(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi,
  435.      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
  436.       qspof=0.0d0
  437.       go to 2000
  438. c
  439. c  finished
  440. c
  441.  1500 qg=0.0d0
  442.       qb=0.0d0
  443.       qc=0.0d0
  444.       qspof=0.0d0
  445.  2000 return
  446.       end
  447.